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Summary 


Two advances in the numerical techniques of utilizing the BIE method 
are presented. The boundary unknowns are represented by parabolas over 
each interval which are integrated in closed form. These integrals are 
listed for easy use. For problems involving crack tip singularities, 
these singularities are included in the boundary integrals so that the 
stress intensity factor becomes just one more unknown in the set of bound- 
ary unknowns thus avoiding the uncertainties of plotting and extrapolating 
techniques. The method is applied to the problems of a notched beam in 
tension and bending, with excellent results. 

1. Introduction 

Knowledge of the stress distribution in the neighborhood of a singu- 
larity, such as the tip of a crack in a beam loaded in tension or bending, 
is of fundamental importance in evaluating the resistance to fracture of 
structural materials. Elastic solutions to various geometries have been 
obtained by a number of different methods. Among the more effective ones, 
are the complex variable method, collocation method, and finite element 
method. However, the first two of these methods are not general enough 
nor readily adaptable to three-dimensional or elastoplastic problems. And 
the finite element method requires solutions of large sets of equations 
and fails to give sufficiently fine resolution in the vicinity of crack 
tips. 

The recently developed boundary integral methods, Mendelson [1], 
offer an attractive alternative to other methods of analysis. These 
methods have a number of advantages which may be listed as follows: 

(1) They obviate the need for conformal mapping. 

(2) Mixed boundary value problems are handled with ease. 

(3) Stresses and displacements are obtained directly without need for 
numerical differentiation. 

(4) No special considerations are needed for multiply connected 
regions. 

(5) The internal stresses and/or displacements are calculated only 
where and when needed. 
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(6) The extension to three-dimensional problems is direct • 

(7) Nodal points are needed only on the boundary instead of through- 
out the interior as required by finite element methods. 

The last point is probably the most important one. For the finite 
element method, the whole region must be covered by a grid producing a 
large number of nodal points and corresponding unknowns. Thus a large 
number of simultaneous equations must be solved. For the boundary inte- 
gral methods, nodal points are taken only on the boundary, resulting in a 
much smaller number of unknowns. On the other hand, the resulting 
matrices are full, whereas in the finite-element methods, for example, the 
matrices are usually sparse and can be more or less banded. 

The method generally used is to divide the boundary C into a number, 
e.g., n, intervals and replace the integrals over C by a sum of n in- 
tegrals over the n boundary interVcils- The unknowns are assumed con- 
stant over each interval. This technique gives good results with rela- 
tively few equations to solve for regions with smooth boundaries. Once 
we introduce geometric singularities, however, such as cracks or notches, 
then a very large number of intervals are required in the vicinity of the 
crack boundary to obtain reasonable accuracy. An improvement can be made 
by assuming linear variations of the unknowns along the boundary inter- 
vals, Riccardella [2]. However, this still results in a relatively large 
number of intervals for reasonable engineering accuracy. 

The present paper presents two major advances in the numerical tech- 
niques of utilizing the HIE method. Firstly, the boundary unknowns are 
represented in terms of parabolas over each interval which are. integrated 
in closed form. These integrals are listed for easy use. Secondly, for 
problems involving crack tip singularities, these singularities are in- 
cluded in the boundary integrals so that the stress intensity factor be- 
comes just one more unknown in the set of boundary unknowns. \-fhen the set 
of linear algebraic equat -ons is solved for the boundary unknowns , the 
stress intensity factor is obtained at the same time, thus avoiding the 
uncertainties of the usual plotting and extrapolating techniques. These 
innovations result in greater accuracy than was possible heretofore, using 
substantially fewer boundary intervals and consequently less computer 
time. 

Results are presented for two practical fracture mechanics configu- 
rations of the edge-cracked plate in pure bending and in tension, and com- 
pared with the standard values quoted in the literature, excellent agree- 
ment being obtained- 

2. Analysis 

Although boundary integral methods can be formulated in many ways, 
for elasticity and elastoplatic problems the most natural formulation is 
in terms of the Navier equations of equilibrium. A solution to these 
equations can be obtained by making use of Kelvin* s singular solution of 
the Navier equations due to a point load and also making use of Betti’s 
reciprocal theorem, Rizzo [3]. We then arrive at a solution that is known 
as Somigliana* s identity, namely: 
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where u^ and Pj are the boundary displacements and boundary loads, re- 
spectively, and tne usual tensur notation is used. The tensors > "^ij* 

are given by 

C3 1 f 

T. . = — ( 6 . . C, + 2r , .r, .) + C, (r, .n. - r, .n.)l 1 

x 3 r j_3ti 13 4 *i ’ 3 ' 4' ’3 x ’i 3 ‘j J 


" 8irG(l - ]i) ’ 


S 4u(l - p) 


= 3 - 4ii 


= 1 - 2u 


and r is the distance from the fixed point P to the variable point of 
integration, s. The above equations are for the case of plane strain. 

For plane stress one replaces Poisson’s ratio y by y/(l + y) . The co- 
efficient X is equal to 1, if P is an interior point and is equal to 
1/2, if P “ p is a boundary point. 

For X equal to 1/2 Eqs, (1) become a set of 2 or 3 (plane problem 
or 3 -d problem) Fredholm equations for the boundary unknowns. These may 
be, boundary tractions, displaceir^nts or combinations of the two. Thus 
the first, second or mixed boundary -value problems of elasticity can be 
solved with equal ease. Once the appropriate unknowns are determined on 
the boundary, the displacement at any interior point P, can be obtained 
from Eqs. (1) with X equal to one. The stresses can be obtained by 
appropriate differentiation under the integral sign of the tensors 
and Tj^j . No numerical differentiation is required. 

The problem then resolves itself the solution of Eqs. (1). The 
method generally used is to divide the boundary C into a number, e, g,, 
n intervals and to replace the integral over the boundary G by a sum of 
n integrals (for each equation) over the n boundary intervals. Tlie un- 
knowns Uj and Pj being assumed constant over each interval. Improved 
results can be obtained by assundng the unknowns to vary linearly over 
each interval, Riccardella [2]. In the present investigation it was as- 
sumed that the unknowns varied parabolically over each boundary interval 
(see appendix). 

When Eqs. (1) are replaced by a sum of integrals over the n inter- 
vals and the integrations carried out, a set of 2 n equations arc 
obtained in 2 n unknowns which can be written in matrix form as 
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where a, 3, a’ , B* are n x n matrices, whose elements are described in 
the appendix, u and v are n x 1 matrices of the x and y boundary 
displacements, respectively, and ^2 U ^ 1 matrices determined 

from the boundary loads as described in the appendix* 

Equation (4) has been written in a form applicable to the first 
boundary value problem of elasticity, i.e*, the tractions are assumed 
known over the complete boundary. For the mixed boundary value problem, 
i.e,, the displacements are given over part of the boundary, the roles of 
those displacements and the corresponding boundary tractions are inter- 
changed in Eq. (4) , This merely involves replacing the appropriate col- 
umns. 

2.2 Stress Intensity Factor 

The usual procedure for determining stress intensity factors for 
specimens with cracks is to calculate the stresses just ahead of the crack 
tip, or the displacements behind the crack tip, or both, and then make use 
of appropriate graphical procedures. For reasonable accuracy it is nec- 
essary, when using numerical methods, to calculate the stresses at a suf- 
ficient Humber of points lying approximately between 0.01 a and 0.1 a 
ahead of the crack tip, where a is the crack length. This will usually 
involve a very fine grid or close interval spacing in the crack tip 
vicinity. 

A method has been developed which avoids both the necessity for a 
large number of intervals in the crack tip vicinity as well as the need 
for graphical procedures with their inherent inaccuracy. The technique 
will be illustrated for the specific problem of an edge-cracked plate 
loaded either in tension or bending (or both) . 

Consider the edge-cracked plate shown in Fig. 1. A distribution of 
end loads producing pure bending is shown, but if desired a tensile load 
could be applied instead. If the complete boundary CDEFAF*E'D^G is 
used, we have a first boundary value problem. However , since in this 
case the problem is symmetric, half the plate can be used if desired, in- 
volving the boundary ACDEFA. We then have a mixed boundary value prob- 
lem, the tractions being unknown on the boundary segment AG, 

Considering the case of the half boundary, we exclude two small in- 
tervals e and in the vicinity of the crack tip. The rest of bound- 

ary is divided into intervals with nodal points taken at their centers as 
shown. Equations (1) and (4) are applied to this boundary with the e 
and intervals excluded. For the intervals e and use is made 

of the known relations defining the mode 1 stress intensity factor, Kj, 
i.e., 



v = 0 


5 


= 4(1 - u") 


P = u = 0 

y 


0 i r ^ e. 


where Py is the traction on the boundary AC, and r is measured along 
the boundary from the orach tip. Substituting into Eqs. (1) and assuming 
Tj j and to be constant at their midpoint values over the Intervals 

ei and e, respectively, gives 
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are evaluated at e/2, and Uj _2 ^.re evaluated at e^/2- 

Adding Eqs* (7) to the rest of the integrals in (1) gives the inte- 
grals over the complete boundary and also adds one more unknown, namely 
Kj, to the 2n unknowns appearing in Eq, (4). Since the number of un- 
knowns is now one greater than the number of equations, an additional 
equation is required. This can be obtained either from the conditions of 
continuity of the boundary tractions at r = e, or the continuity of the 
crack opening v at r = or both. Using the parabolic assumption 
gives the equations 
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It was found that the best results were obtained by either using both 
equations giving a slightly overdetermined system, or by adding tha two 
equations together to give a single equation. 


3. Results and Conclusions 


Calculations were carried out for the plate shown in Fig* 1 for both 
bending and tension. The relative dimensions of the plate were taken as 
W = 1, a= 0.5, L= 1.2, with p = 0.33. In all the calculations each of 
the segments CD, DE, and EF was divided into seven intervals. The numbers 
of intervals along FA and AC (NFA and NAG) were varied to determine the 
effect of increasing the number of these intervals. 


Some of the results are shown in Figs. 2 to 5. The crack opening 
displacements and the nornial stresses ahead of the crack are plotted and 
compared with results obtained by Gross [4] using boundary collocation. 
Excellent agreement is obtained using relatively few nodal points. As a 
matter of fact taking NAC = NFA = 10, changes the values to 1.76 and 

3.46 for the cases of bending and tension, respectively. This is still 
sufficiently accurate from an engineering viewpoint and makes use of a 
total of 41 boundary intervals resulting in just 83 unknowns (including 
Ki) . 


It can be concluded from the results obtained that the Boundary Inte- 
gral Equation method utilizing the techniques described herein can be an 
efficient and accurate tool in fracture mechanics analysis. 

4. Appendix 


4,1 Parabolic Assumption 

It is assumed that the unknowns vary parabolically over each boundary 
interval. Thus for the jth boundary interva:! 
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where Uj and Vj are the displacements in the x and y directions, 
respectively, and s is measured from the beginning of the interval under 
consideration. Similar formulas apply to the boundary tractions. 


The coefficients in Eq, (10) are given by 
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with similar relations for hQj , bj^j, and b2j* h is the interval length 

which in the above formulas, is assumed to be the same for the jth inter- 
val and the adjacent 3 - 1 and j + 1 intervals. Similar type formulas 


can be obtained If desired for unequal intervals. 


If the interval under consideration is an end interval such as, for 
example, the first interval of segment AC and the last interval of the 
segment CD of Fig. 1, then Eqs. (11) are modified as follows. For a begin-- 
ning interval such as interval 1, 
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For an end interval such as interval ND 
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4.2 llatrix Coefficients 

The coefficients entering into the matrices of Eq. (4) are obtained 
as follows. The boundary is divided into n intervals and Eqs. (1) 
replaced by sums of integrals over the boundary intervals. Equations (2) 
and (11) through (13) are substituted in and the integrals evaluated in 
closed form. Let i designate a fixed nodal point on the boundary and j 
a boundary interval whose length is ASj as shown in Fig. 6. r j , ^j*fl> 

0j , and D are defined as shown in the figure. Note that 

D = r j cos 0j = rj^.^ cos Using this notation we define a set of E 

matrices whose elements of the ith row are given by 
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Also let JZ.. and m. be the direction cosines of the normal n. to the 
J J J 

interval, and define 
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Then the following matrix coefficients are computed. 
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Note that for the problems under consideration, the last six integrals in- 
volve intervals on side AC for which X, = 0 and m = -1. 


The vectors and P 2 appearing in Eq. (4) are given by 
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Carrying out the integrations gives 
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where: T = sin 
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For tension: 
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Tlie matrices a, a' , 8, 8' of Eq. (4) are now assembled from the 
matrices A through K' using the parabolic coefficients of Eqs . (11) to (13) 
As an example, consider the line AC having NC nodal points with interval 
length HAC. Then for 3 < j < NC - 3 , the central formulas are used and 
we write 
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At the end of the line we have 
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Exactly similar relations hold for a’, P, and S’, replacing A, G, and 
D by (A’ , C’ , D’), (B, E, F; and (B’ , E’ , F'), respectively. 
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